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INTRODUCTION 

Transient thermography employs preheated surface temperature variations caused by 
delaminations, cracks, voids, corroded regions, etc. Often, it is enough to detect these 
changes to declare a defect in a workpiece. It is also desirable to obtain additional 
information about the defect from the thermal response. The planar size, depth, and thermal 
resistance of the detected defects are the parameters of interest. 

In this paper a digital image processing technique is applied to simulated thermal 
responses in order to obtain the geometry of the inclusion- type defects in a flat panel. A 
three-dimensional finite difference model in Cartesian coordinates is used for the numerical 
simulations. Typical physical properties of polymer graphite composites are assumed. Using 
different informative parameters of the thermal response for depth estimation is discussed. 

NUMERICAL MODEL DESCRIPTION 

A wide variety of the one- two- and three-dimensional models have been applied to 
investigate the capabilities of thermal NDE for defect characterization. The parabolic heat 
diffusion equation provides a broad basis for analytical and numerical modeling of the heat 
propagation in a solid body. The one-dimensional approach allows study of the transient 
thermography task, with the assumption that the lateral size of a subsurface defect is big 
compared to the thickness of the material above it. The two-dimensional model offers more 
flexibility for the defect configuration. Using this model, the thermal response from a defect 
with finite width and thickness can be obtained and studied. It is also possible to monitor a 
cross influence of two or more defects. Unfortunately, only one-dimensional surface 
temperature profiles can be produced by a two-dimensional model. This limits research of the 
image processing techniques which can be applied to the thermal images obtained from an 
infrared camera. 



2 


To estimate capabilities of the thermal method for quantitative characterization of 
inclusions, a numerical model of the three-dimensional heat transfer problem has been 
developed [1]. The influence of the values of the different defect parameters (planar shape, 
depth, and thermal resistance) on the thermal response of a composite panel is investigated. 
The finite difference mesh has elements of size 1 mm in the xy -plane and 0.2 mm in the z- 
direction. It represents a rectangular, thermally insulated, plate with one, two, or more thin 
rectangular inclusions of variable size. A typical geometry of the defects in the plate is shown 
in Fig. 1. The energy balance for the control volume is used on the defect-plate boundaries. 
For carbon-epoxy composite plate the thermal diffusivity value of 2.8 x 10” 6 m 2 / 5 is used. 
Thermal properties of the inclusions are varied to represent the defects of different nature. 

Pulse heating of the upper surface with coordinate D_ is used. The heat is applied 
during only the first time step and can be considered as instantaneous. One-sided access to 
the top surface is considered. The finite difference program uses the simple explicit method 
to solve the three-dimensional transient heat exchange problem. The result of the program is 
the temperature of each node of the mesh at each time step. 

It is easy to produce a surface thermal response as well as extract a time evolution 
curve for a node from the simulated data. Designed in such way, the program can be easily 
reduced to the two- and one-dimensional problem by expanding the defect to the plate edges. 
In addition, it has an ability to produce the thermal response from a group of two or more 
defects. This provides a basis for comparing the simulated data and the infrared thermograms 
obtained from a physical model of the defects. 

INVERSE TASK 

When the physical properties of the inclusions are equal to the properties of the 
surrounding material, there are no spatial variations of temperature on the top surface of Fig. 
1. The heat flux becomes one-dimensional, and the model behaves as an opaque uniform 
sheet of thickness D z . The physical parameters of the inclusions and their position cause 
changes of temperature on the surface compared to the unflawed state. An inversion 
procedure characterizes the unknown defects from the temperature evolution on the plate 
surface after heating. Reconstruction of the planar shape of the delaminations and their depth 
from the received temperature patterns is considered in the following sections. 



Figure 1. Model of a plate with the defects for the numerical simulations. 



Planar Size Extraction 
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Different methods of surface thermal response processing can be used to determine 
the planar size of a defect. One of them is a calculation of the spatial gradient of the surface 
thermal contrast. The distance between the steepest parts of the contrast above a defect has 
been reported as a stable parameter for the defect size estimation [2]. This approach makes 
unnecessary an earlier estimation of the threshold level and can be applied to a thermogram 
which contains several defects. The second partial derivative of the contrast C is proposed, 
as explained later, to compute this distance. The contrast between the temperature above a 
defect T def and temperature above the defect free (sound) area T soa is calculated as 


C(t) = 


^def (0 T soa (t) 


T soa { t) 


( 1 ) 


A model with two inclusions of different depth and size is used to produce a 
temperature pattern on the surface. It allows the simulation of a general class of specimens 
with unknown number, size and location of the defects in the observed area. The inclusions 
are square air defects (thermal resistance R — 2.9 x 10 _3 m 2 AT/ W ) with a size of 10 and 20 
mm on a side, located 2.4 and 7.4 mm respectively under the surface of the 10 mm-thick 
plate (Fig. 1). 

To obtain the thermal response from a defect located at any depth, the contrast is 
averaged during the observation time interval of 36 s. The length of the interval is chosen 
from the condition when the Fourier number (normalized time) equals one. The 3D contrast 
distribution for the defects is shown in Fig. 2a. Two peaks of magnitude above the defects 
are noticeable in the plot, as well as the cross influence of the defects. The contrast profile 
across the center line of the defects is shown in Fig. 2b. The profiles of the averaged contrast 
for Teflon inclusion (R = 0. 11x10 ~ 3 m~K / W ) and the low conductive defect 
(R = 200 x 10 3 m 2 K !W) are also presented. As can be seen from this figure, the steepest 
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Figure 2. The contrast above two square defects, a) three dimensional plot of the contrast for 
air inclusions; b) thermal contrast curve for different thermal resistance of the inclusions. 
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Figure 3. Planar defect size extraction using the second derivative method, a) Illustration of 
the defect detection algorithm; b) directions for the algorithm. 


parts of the profiles keep approximately the same location for various depth and thermal 
resistance. 


The first and second x-derivatives of the contrast profile above the air inclusions are 
sketched in Fig. 3a. The parts of the contrast curve where 


dx 2 


<0 


( 2 ) 


define the size of the defects close to their modeled size (shaded in Fig. 3a). An application 
of this criteria in the opposite direction provides the same size of the defects. To reconstruct 
an 



e) f) g) h) 

Figure 4. Defect maps, a) i-map; b) j-map; c) i+j+ diagonal map; d) i+j- diagonal map; e) 
result of multiplication of the maps (a) and (b); f) result of multiplication of the maps (c) and 
(d); g) sum of the maps (e) and (f); h) thresholded image of (g). 
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unknown shape of the defects on the rectangular grid, a finite difference approximation of 
this partial derivative is used. Suppose, that increment of i and j represents the positive 
direction of axes Ox and Oy respectively. To make a decision about defectiveness in the 
pixel with the coordinates (i,j) the 8 surrounding nodes are considered (Fig. 3b). The criteria 
(2) can be applied in the perpendicular or diagonal directions giving “1” for a node which 
satisfies the criteria and “0” if it does not. 

The diagrams in Fig. 4 illustrate the application of this procedure to the summed 
contrast image shown in Fig. 2a. The defect map computed using the criteria (2) for the 
nodes (i - 1 , j ) , (/,/), and (i + 1 , j ) (i-map) is shown in Fig. 4a. The map identifies both 
defects, but 

gives incorrect size in the j -direction. Similarly, the j -map (Fig. 4b) produces incorrect shape 
of the defects in the i-direction. Result of multiplication of these two binary maps (Fig. 4e) 
provides the correct reconstruction of the planar shape of both square defects. Because a 
defect image on an experimental thermogram can have a non-orthorhombic configuration, it 
seems reasonable to consider also the diagonal directions for the defect contour extraction. 
The defect maps received by the application of this procedure at the diagonal directions 
(i+,j+) and (i+,j-) (Fig. 3b) are shown in Fig. 4c-d. The defects on the image derived from 
these maps have rounded corners (Fig. 4f). The next step of the defect imaging process is a 
combination of the orthogonal and diagonal maps. The sum of these maps is shown in Fig. 
4g. The thresholded defect image is in Fig. 4h. This alternating direction algorithm is 
expected to be more robust in practical applications, although for the our model task it 
produces worse (more rounded) reconstruction of the defect shape. 

An advantage of the proposed algorithm is that it can be applied to a segment of a 
thermogram. To illustrate this statement, the contrast image was divided into three segments 
(Fig. 5a). Segment 1 covers only a quarter of the second (deepest) defect. Segment 2 contains 
one quarter of the second defect in the region where its image is affected by the first defect. 
Segment 3 crosses both defects at the centerline. The algorithm was applied independently to 
each image segment. The three computed binary images were linked into one defect map 



a) b) 

Figure 5. Application of the defect map construction algorithm to three segments of the 
simulated thermogram, a) gray scale image of the thermal contrast above two defects (see 
text); b) defect map received after a linkage of the segments. 
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(Fig. 5b). The solid-line contours of the defects are also presented for comparison. The image 
is different from the Fig. 4h even though it keeps the position of the defects and their 
respective geometry. One concern of this technique is the case of a defect of large lateral 
dimension such that the method will produce only a wide contour of the defect since the top 
is flat and the spatial derivative is zero. 

Depth Estimation 

Several informative parameters of the modeled thermal response have been calculated 
to study the depth estimation task. They include the maximum value of the thermal contrast 
Cmax and the time t max when this maximum occurs, the time of the peak slope t of the 

thermal contrast curve, the amplitude and the phase of the Discrete Fourier transform applied 
to the temperature evolution on the surface after the heat injection. These informative 
parameters are calculated at the node which belongs to the surface and is located above the 
center of the subsurface defect. The size of the plate used for simulations is 80 mm x 80 mm 
x 20 mm. The defect is 0.2 mm thick air inclusion with the planar size of 20 mm x 20 mm 
(normalized size L = 1/ D_ = 1). It is located at the center of the plate. The depth of the defect 
under the surface varies from 1 mm to 19 mm. The calculated parameters are shown in Fig. 

6. As seen from these plots, all informative parameters are sensitive to the depth variation. 

The amplitude and phase values are calculated from the imaginary and real 
components of the result of the Discrete Fourier transform computed for the fundamental 
frequency at a time interval of 75 s. They have non-homogeneous character and are not 
applicable for the defect depth extraction. The three other parameters have been actively 
studied for depth estimation both theoretically and experimentally. It has been found that the 
temporal characteristics t max and t are less sensitive to the defect size and thermal resistance 

variations than C max . 

Several authors report, that the time of the peak slope is a steady characteristic of the 
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Figure 6. A comparison of the informative parameters for the air defect with size of 20 mm x 
20 mm x 0.2 mm. 
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Figure 7. Time of the peak slope t for various defect size and thermal resistance. 

defect depth. For example, the experimental research of this parameter has been reported by 
Favro et al. [3]. The experiments were accomplished with the flat bottom holes of different 
diameter and depth. The t magnitude, as reported, rises as a square of the defect depth and 

has a low sensitivity to hole diameter variations. The inverse task of depth estimation can be 
solved by taking the square root of the time, measured experimentally. In contrast, the 
t ps curve plotted in Fig. 6 has a relatively linear increase. A detailed analysis has been 

performed to find an agreement between these data. Fig. 7 contains a comparison between 
the results received on a one-dimensional model and on the 3D finite difference model. For 
the one-dimensional case two infinite plates of different thickness are used. One represents 
the material above a non-conductive defect. The second, with thickness D_ , is used as a 
reference. The analytical solution of the heat diffusion equation for the surface temperature 
of an insulated plate after a Dirac heat pulse [4] is used. The thermal contrast is calculated as 
a difference between the temperatures on the plate with the variable thickness (representing 
the defect depth) and on the plate with the reference thickness. 

Three t curves received on the one-dimensional model for D : of 20 mm, 40 mm, 

and 100 mm are shown in Fig. 7. When the thickness of the reference plate is big enough 
compared to the depth range, the t function has a parabolic shape (curves 2 and 3 in Fig. 7). 

The t ps curve received for D_ = 20 mm diverges from the parabola at the middle of the depth 
range. All values of t received from 3D modeling fall below this one-dimensional case 

result. The time increases as the defect size or thermal resistance increases. For instance, the 
curve for the low conductive defect (R = 200 x 10 ~ 3 m 2 K /W ) with normalized defect size 
k = 2 lays closely to the one-dimensional result. It is seen that the planar size and thermal 
resistance of the defect both have a strong influence on the peak slope time, which decreases 
accuracy of the depth estimation based on this parameter. 


CONCLUSION 
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A convenient algorithm for defect map construction using the distance between the 
steepest parts of the surface thermal contrast above a defect for the size calculation is 
presented. This method can be applied to a thermal image without preliminary localization of 
the defects. The robustness of the method has been tested on image segments containing only 
parts of the defects. The described defect mapping procedure will be verified on a physical 
model of the defects. 

By way of numerical simulations a complex dependency of the time of the peak slope 
of the thermal contrast on the defect planar size and thermal resistance has been discovered. 
As followed from the received results, the usage of the square root of this time for defect 
depth calculation is limited to the depths less than a half of the plate thickness. 
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